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ABSTRACT. In serine proteases (SP's), the H-bond between 
His-57 and Asp-102, and that between Gly-193 and the transition 
state intermediate play a crucial role for enzymatic function. To 
shed light on the nature of these interactions, we have carried out 
ab initio molecular dynamics simulations on complexes represent- 
ing adducts between the reaction intermediate and elastase (one 
protein belonging to the SP family). Our calculations indicate the 
presence of a low— barrier H-bond between His-57 and Asp-102, in 
complete agreement with NMR experiments on enzyme— transition 
state analog complexes Comparison with an ab initio molecular 
dynamics simulation on a model of the substrate— enzyme adduct 
indicates that the Gly-193— induced strong stabilization of the in- 
termediate is accomplished by charge/dipole interactions and not 
by H-bonding as previously suggested. Inclusion of the protein 
electric field in the calculations does not aff'ect significantly the 
charge distribution. 

KEY WORDS: density functional theory calculations; enzyme- 
intermediate adduct; H-bonding interactions; low— barrier hydro- 
gen bonds; Car— Parrinello simulations 
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INTRODUCTION 

Hydrogen bonding is one of the most important interactions for the bi- 
ological function of enzymes. Extensive H-bond networks may stabilize the 
conformation of the active site which is necessary for the catalytic function 
and they can fix the optimal orientation of the substrate for the enzymatic 
reaction. Most importantly, H-bonds may allow for transition state stabi- 
lization by lowering the activation free energy by several kcal/mol [0, ^. 

A typical example in this respect is constituted by serine proteases (SP's) 
enzyme family §, |, §, §, 0, |, ||, |TU|, [TT], 0. 



SP's use the catalytic triad (Ser-195-His-57-Asp-102) to catalyze the hy- 
drolysis of peptides (Fig. la). This occurs through nucleophilic addition of 
the 3-hydroxyl group of Ser-195 to the acyl carbonyl of the substrate, with 
formation of a negatively charged tetrahedral intermediate (Fig. lb). 

Stabilization of the intermediate is achieved by formation of two H-bond 
with the amide groups of Ser-195 and Gly-193 (mammalian isoenzymes |^]) 
or with the amide groups of Ser-195 and the sidechain of Asn-155 (bacterial 



iso enzymes |Tj] ) 



Theoretical [T^, |T6| and experimental [0, [T^ studies on wild type and 



mutants of a bacterial SP (subtilisin) have shown that Asn-155 is a key 
residue for the biological function, in that it provides a stabilization of the 
transition state (TS) relative to the ground state (GS) by as much as ~ 5 
kcal/mol. Curiously, no correspondent studies on the mammalian isoenzymes 



have appeared to shed light on the crucial role of Gly-193|]l8 . 

A second, important H-bond interaction involve two residues of the cat- 
alytic triad, His-57 and Asp-102. A series of NMR studies on a mammalian|l|. 



Til [H) mi bacterial ||2^ SP's and their complexes with inhibitors have 
indicated the presence of a low-barrier hydrogen bond (LBHB) linking N^i 
of protonated His-57 with the P carboxyl group of Asp-102 (Fig. lb)|l|, 

[T^. Approaching of the TS is suggested to facilitate the formation of the 
LBHB, which in turn may render N^2 of His-57 a stronger base for accepting 
a proton from Ser-195 in the formation of the intermediate [0, |TT], |T2|, |TB . 
As result of this process, the free energy barrier of TS relative to GS de- 
creases. Ab initio calculations and neutron scattering experiments have led 
to the conclusion that the LBHB is very covalent in nature PT|. Thus, this 



"resonance-stabilization" energy could supply much of the energy necessary 
for enzyme's catalysis |]T], [IT], |T3|. However, the role of this LBHB for the 
catalytic power of SP's is object of controversy P^, |25[ . 



From the above considerations it is clear that, in spite of its crucial role, 
the nature and the dynamics of hydrogen bonding in active site of SP are not 
fully understood. In order to provide a picture of the chemical bonding of 
these interactions in SP's, and to relate it to the biological function, we have 
carried out ab initio molecular dynamics simulations on models of SP- 



intermediate and SP-substrate complexes. This technique is revealing itself 
as a useful computational tool for investigating specific molecular interactions 
in biological systems [^]. Based on state-of-the art density functional theory 
calculations, it describes interactions and charge distributions of (relatively) 
large systems rather accurately. Furthermore, it includes temperature effects, 
which obviously play a fundamental role in SP function. Finally, they allow 
describing bond breaking and bond forming processes at room temperature, 
which is essential to describe phenomena such as the LBHB. 



Anticipating our results, the calculations show that the LBHB is a strongly 
covalent interaction and that the dramatic Gly-193-induced stabilization of 
the reaction intermediate is due mainly to the electrostatic interactions be- 
tween the intermediate and the Asn-192-Gly-193 peptide's unit dipole. 



COMPUTATIONAL PROCEDURE 
Model complexes. 

Our structural models for the adducts of SP with intermediate (I-SP) and 
substrate (S-SP) are based on the X-ray structure of porcine pancreatic ela- 
stase (240 residues) complexed with Ace-Ala-Pro- Val-difluoro-N-phenylethylacetamide 
(PDB entry: 4EST) |26|, . They include the entire catalytic triad, the scis- 
sile peptide bond and the oxyanion hole. 

The construction of the complexes is carried out in several steps: i) the 
terminal N-phenylethylacetamide is replaced by an acetyl group; ii) all hy- 
drogen atoms of the complex, not present in the X-ray structure, are added 
assuming standard bond lengths and bond angles; in) a shell of 1453 water 
molecules including the crystallographic ones (approximatively correspond- 
ing to three water monolayers), is added; iv) four chlorine counter ions are 
added to ensure neutrality; v) energy minimization is carried out with the 
AMBER suite of programs using the AMBER force field |2|] (conver- 
gence criterion 0.0001 kcal/(mol-A)). In the minimization, no periodic bound- 
ary condition are applied and the electrostatic interactions are calculated 
assuming a constant dielectric function e = 1 and without cutoff. 
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The I-SP and S-SP model complexes comprises the entire side-chain of 
Asp-102, the imidazole ring of His-57, the Q192G193 peptide linkage and 
the entire Ser-195 residue, which in I-SP is covalently bound to the substrate 
(Fig. ^,b). His-57 is considered doubly protonated in I-SP and protonated 
in the 6 position in S-SP . 



Alternative models representing I-SP and S-SP differ from the previous 
ones for the substitution of the Q192G193 peptide link with dimethylammo- 
nia (Fig. 0c, d). 

Calculations. 

The quantum-mechanical problem is solved within the framework of den- 
sity functional theory (DFT)[Q, in its Kohn and Sham (KS)^1|] formulation. 
The KS orbitals are expanded in a plane- wave basis set, up to an energy 
cut-off of 70 Ry. Only valence electrons are considered explicitly, pseudopo- 



tentials of the Martins and Troullier p2| being used for the core-valence 



electron interaction. BLYP p3|, [3^ gradient-corrected exchange-correlation 



functionals are used. The charge of all the complexes is -1. 

We have carried out geometry optimization using the direct inversion in 
iterative subspace (DIIS) method ^ for both electronic (convergence 
threshold set to 10^^) and ionic degrees of freedom (convergence threshold 
set to 5 - 10-^). 



DFT-based Car-Parrinello ab initio molecular dynamics simulations 
are performed at constant temperature, with the atomic trajectories col- 
lected over a period of 0.6 and 1.0 ps for S-SP and I-SP, respectively. Equa- 
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tions of motion are integrated with the velocity Verlet algorithm. The ficti- 
tious electron mass is 400 a.u. and the integration timestep is 4 a.u. Con- 
stant temperature simulations are achieved by coupling the systems to a 
Nose thermostat [^] at T = 300 K with a frequency of 500 cm~^. The ter- 
minal hydrogens of Asp-102, His-57, and Gly-193, corresponding to the C^, 
Cjs and 0^ respectively, are kept fixed in their starting position; in S-SP an 
additional constraint between the 0^(Ser-195) and the substrate carbon of 
the scissile bond (indicated in Fig. |^ as Cg) is imposed. 

Calculations including the external electrostatic potential of the whole 
protein-water system are also carried out. This potential ^protir) at the 
point r is evaluated as 

z — ^ \r- — r\ 



where qi are the RESP I^H] atomic point charges at point r,. 

All Car-Parrinello calculations are performed with a parallel version of 



the CPMD code V3.0h i39f 



Calculated properties. 

The electrostatic energy AE between to moieties (e.g. the intermediate 
and Q192G193 peptide unit) is calculated as 

AE = J2^-^, (2) 



where the indexes i and j refer to atoms of the two moieties, qi and qj are 
the partial atomic ESP charges PD[| and Vij is the interatomic distance. 



Test calculations are carried out also using the multipolar expansion of 
the electrostatic energy up to the dipolar term: 



AS ~ Qi<l>2 - Ml ■ E2. (3) 

where Qi and ^tl are charge and dipole moments of moiety 1 and $2 
and E2 the electric potential and the electric field produced by moiety 2. 
respectively. The results turn out to be very similar to those obtained with 
eq. i 

Binding energies (B.E.'s) are calculated as total energies differences be- 
tween complexes in Fig. ^and their forming elements. The B.E. of complexes 
I and III could not be determined because of the instability of the interme- 
diate fragment. 



RESULTS 

In this section, first we analyze structural and electronic features of two 
models representing the adduct between serine protease and the reaction 
intermediate (I and III of Fig. Comparison is then made with features 
of models of the substrate-enzyme complex (II and IV of Fig. 

The Intermediate— Enzyme complex 

Structural features and charge distribution. Conformational properties 
as well as the H-bond network of the complex are fairly maintained dur- 
ing the dynamics (Fig. |^). Consequently, the charge distribution does not 
change significantly (with exception of C5 and N5), as indicated by the ESP 
atomic partial charges reported for several snapshots of the molecular dy- 
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namics (Tab. 1). Note that the C5-O bond of the intermediate is very 
polarized towards the oxygen, consistently with the fact that this bond is 
to be broken in the subsequent step of the hydrolysis. The presence of the 
protein field does not affect significantly the charge distribution (Tab. 2), 
suggesting that solvent effects do not play a major role for the electrostatic 
interaction at the active site. 

H-bond pattern: Asp-102-His-57. During the dynamics, proton hopping 
occurs between the His-57 and Asp-102 in the subps time-scale (Fig. 
The presence of a LBHB is completely consistent with NMR data on an 
intermediate-serine protease complex, namely the peptidyl trifluoromethyl 
ketone-chymotrypsin adduct []T|. 

The chemical bonding of the LBHB can be characterized with the electron 
localization function (ELF) ^ The ELF has proven to be very useful 



to illustrate chemical concepts like localized bonds and electron lone pairs. 
Fig. ^ shows the ELF before, during and after the proton transfer from His- 
57 to Asp-102. The red areas indicate strong localization, i.e. spatial regions 
where the Pauli principle has little influence on the electron distribution. 

Fig. I^a shows the presence of the lone pairs of the aspartate oxygen and 
of the strong electron localization along the histidine N^i-H bond, which 
indicates the covalent nature of the bond. During the proton transfer (Fig. 
|b), the N^i-H bond is still very covalent and an incipient covalent 0^2^13 
bond is being formed. Protonation of Asp-102 establishes a covalent 0-H 
bond: significant portion of ELF is indeed localized on the H atom (Fig. §c). 
The formation of the nitrogen electron lone-pair is also evident from the 
picture. We conclude that the bonding in this LBHB is essentially covalent 
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in nature. Similar findings have been reported in a very recent ab initio study 
of low-barrier H-bonds in an organic molecule ||21|| . 

H-bond pattern: Gly- 1 93-intermediate. The second fundamental H-bond 
interaction investigated here involves Gly-193 and the intermediate carbonyl 
oxygen. This H-bond is well maintained during the dynamics (average O- ■ -H 
distance of 1.7(0.1) A). A rough estimation of the interaction energy based 
on an electrostatic model (see Computational Section) indicate that Gly-193 
stabilizes the intermediate by more than 10 kcal/mol (Tab. 3). This value 
appears to be too large for a purely electrostatic H-bond interaction ^ 



Inspection of the structure reveals that the very large Q192G193 peptide's 
unit dipole (;^ 4 D [^) could be also an important factor for intermediate 
stabilization, as it points towards the negative charge of the intermediate. To 
extract the peptide dipolar contribution from the total stabilization energy 
we construct a second model complex in which the Q192G193 peptide unit 
is substituted by dimethylammonia (III in Fig. ^). Tab. 3 shows that 
the resulting stabilization is much smaller, only few kcal/mol. Thus, we 
conclude that a large contribution of the transition state stabilization is due 
to electrostatic interaction (charge-dipole interactions). 

To study the relevance of the Q192G193 dipole on the dynamics, an ab 
initio molecular dynamics simulation on complex III, where the Q192G193 
peptide unit is replaced by a dimethylammonia, is performed. Fig. ^ which 
reports structural properties of the complex, indicate that this complex is 
very unstable with respect to the substrate-enzyme complex. Indeed, while 
the key Gly-193-intermediate H-bond becomes very weak(Fig. |]a), the proto- 
nated His-57 transfers a proton to the intermediate (Fig. ^a) and the 0^(Ser- 
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195)-C5' bond of the intermediate breaks. As a result, a double 05-0(1) bond 
is formed (as indicated by the decrease of the bond distance up to the typi- 
cal value of a carbonyl peptide bond (1.25 A in Fig. ^b)) and 0^ changes it 
hybridization from sp^ to sp^, with formation of the planar peptide unit (as 
shown by the increase of the Z(N(I)-05-0(I)) angle up to ^ 120° (Fig. ||c)). 
In conclusion, our calculations suggest that the absence of the stabilizing 
Q192G193 dipole causes the reverse of the reaction, with formation of the 
substrate and the original H-bond pattern of the catalytic triad. 

The Substrate— enzyme complex 

Gin- 1 92-substrate interactions. To estimate the stabilization of the Q1920193 
peptide unit's dipole on the substrate, we perform an ab initio molecular dy- 
namics simulation of a model of the enzyme-substrate adduct (II in Fig. ^). 

Fig. 7 shows that during the dynamics the two key H-bond interactions 
are maintained but no proton transfer occurs. Interestingly, the substrate- 
protein interaction energy turns out to be much lower than that of the I-SP 
complex (Tab. 3). Replacing the Q192019 peptide with dimethylammonia 
(complex IV) causes a drastic decrease of the interaction energy. The latter 
turns out to be practically identical to that of complex III (Tab. 3). We 
conclude that the H-bond interaction are similar in the S-SP and I-SP com- 
plexes. In contrast, the electrostatic (charge-dipole) interactions are very 
different, the I-SP being more stable by ~ 6 kcal/mol than S-SP (Tab. 3). 
It is interesting to note that this value compares well with previous quan- 
tum mechanical calculations for the Asn-155-TS stabilization in the bacterial 
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isoenzyme 



T5|, We conclude that the transition state stabihzation is due 



mostly to charge-dipole interactions. 

For these complexes it has been possible to calculate also the binding en- 
ergies. Tab. 3 shows a qualitative agreement between binding energies and 
energies based on electrostatic model. This validates the use of the electro- 
static model for a qualitative analysis of intermolecular interactions, as it has 
been done in this work. However, it must be stressed that use of more real- 
istic quantum mechanical models, which comprises other aminoacids present 
in the active site cavity, is expected to screen the charges and therefore to 
reduce the calculated binding energies. 

Charge Distribution. Also in this complex, most of the ESP charges do 
not vary significantly during the dynamics and by introducing the electric 
field of the protein (Tabs. 4 and 5). Most of the ESP charges turn out to be 
similar to those of the I-SP complex. A notable exception is represented by 
the C-0 peptide bond, which in this case is much less polarized toward the 
oxygen. Thus it appears that the protein active site, and in particular the 
Q192G193 moiety, is engineered so as to render the scissile bond more polar 
in the formation of the transition state. 



DISCUSSION 



Within the very short time-scale here explored, our ab initio molecular 
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dynamics simulations help elucidate important aspects of two key interac- 
tions in serine proteases-reaction intermediate complexes, the His-57-Asp- 
102 LBHB and the Gly-193-intermediate H-bonds. 

Our calculations are completely consistent with and confirm the exis- 
tence of a LBHB between His-57 and Asp-102, which has been observed 
experimentally in transition state analog inhibitor complexes [12], |TB . 



Furthermore, they strongly support the proposal of an LBHB-facilitated 
mechanism |T]]. Indeed, the LBHB turns out to be mostly covalent in nature. 
The energy supplied by covalent interaction may be crucial to overcome the 
energy loss due to the compression of the two residues, which is a prerequisite 
of the postulated LBHB-based reaction 

The second conclusion of this paper is that the rather large, Gly-193- 
induced stabilization of the transition state with respect to the fundamental 
state is not caused by an H-bond with Gly-193, as commonly proposed|Q, 
indeed, as the H-bond favors the binding of both substrate and intermediate 
by ~ 2.6 kcal/mol, a value typical of a strong H-bonds in biological sys- 
tems Instead, the negatively charged transition state turns out to be 



more stable relative to S-SP by several kcal/mol as a result of the interaction 
of the negative charge with the large dipole of the Q192G193 peptide unit. A 
simulation in which dimethylammonia replaces the Q192G193 peptide unit 
confirms the crucial role of the dipole: the absence of the stabilizing charge- 
dipole interaction renders the intermediate species unstable. These consider- 
ations suggest that site-directed mutagenesis experiments on the 192 and/or 
193 positions might affect significantly the activity of SP's, as the Q192G193 
dipole orientation may be not optimal for transition state stabilization. 
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Because environment effects may be very important for the chemistry 



of the active site of this and other enzymes |T5|, [16|, we carry out 

some of the calculations in the presence of the electric field of the protein. 
Our results, summarized by tables 1-2 and 4-5, indicate however that the 
field appears not to affect dramatically the charge distribution of the I-SP 
and S-SP complexes. More sophisticated models of the protein electric field, 
which for instance include the electronic polarizabilities of the protein atoms, 
are not expected to alter significantly the picture. 
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Figure 1: Schematic views of the H-bond network in mammahan serine 
proteases active site (a) and of the adduct with the intermediate of the enzy- 
matic reaction (b) . In (b) the double arrow symbol refers to the a low-barrier 
H-bond. 

Figure 2: Model complexes representing I-SP ((a) and (c)), S-SP ((b) and 
(d)). In (c) and (d) the Q192G193 peptide unit is replaced by dimethylam- 
monia. H-bonds are depicted with dashed lines. Green arrows indicate the 
scissile carbon atom C5. The latter is labeled only in (b) for clarity. 

Figure 3: Molecular dynamics of I-SP: final structure of model I. H-bonds 
are represented with dashed lines. 

Figure 4: His-57-Asp-102 H-bond in I-SP (complex I): H-Oj2(Asp-102) (red 
line) and H-N5i(His-57) (blue line) distances plotted as a function of time. 

Figure 5: His-57-Asp-102 proton transfer: electron locahzation function 
(ELF) of three snapshot during the dynamics. The ELF is represented in a 
best-fit plane containing the oxygen, the proton and the imidazole ring, and 
it ranges from (blue) to 1 (red). 
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Figure 6: Molecular dynamics of I-SP: selected structural properties of Com- 
plex III plotted as a function of time, (a) H(Gly-193)-0(I) (blue line), 
He:2(His-57)-0^(Ser-195) (green line) distances; (b) C5-0(I) bond length; 
(c) N(I)-C5-0(I) angle. 

Figure 7: H-bonding of S-SP (Complex II): H(Gly-193)-0(S) (blue line), 
O52(Asp-102)-H5i(His-57) (red line) distances plotted as a function of time. 
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TABLE 1. I-SP in vacuo 
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Table 1: Selected ESP partial atomic charges of I-SP in vacuo. 



TABLE 2. I SP with protein external field 
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Table 2: Selected ESP atomic charges of I-SP in the presence of the protein electrostatic potential. 



TABLE 3. Intermediate - and substrate - Q192G193 peptide unit interactions 
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Table 3: Energies (kcal/mol) are calculated from the electrostatic ESP-based model and from binding energies (see 
Computational Section). ESP-based energies of complexes I and II are calculated as average during the dynamics, 
whereas those of complexes III and IV from the initial structural model. 



TABLE 4. S-SP in vacuo 
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Table 4: Selected ESP atomic charges of S-SP in vacuo. 



TABLE 5. S SP with protein external field 
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Table 5: Selected ESP atomic charges of S-SP in the presence of the protein electrostatic potential. 
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Fig. 1. De Santis L., Carloni P. 
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Fig. 4. De Santis L., Carloni P. 
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Fig. 7. De Santis L., Carloni P. 



